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Two mechanisms for nonlinear mode saturation of the r-mode in neutron stars have been sug- 
gested: the parametric instability mechanism involving a small number of modes and the formation 
of a nearly continuous Kolmogorov-type cascade. Using a network of oscillators constructed from 
the eigenmodes of a perfect fluid incompressible star, we investigate the transition between the two 
regimes numerically. Our network includes the 4995 inertial modes up to n < 30 with 146,998 
direct couplings to the r-mode and 1,306,999 couplings with detuning< 0.002 (out of a total of 
approximately 10 9 possible couplings). The lowest parametric instability thresholds for a range 
of temperatures are calculated and it is found that the r-mode becomes unstable to modes with 
13 < n < 15. In the undriven, undamped, Hamiltonian version of the network the rate to achieve 
equipartition is found to be amplitude dependent, reminiscent of the Fermi-Pasta-Ulam problem. 
More realistic models driven unstable by gravitational radiation and damped by shear viscosity 
are explored next. A range of damping rates, corresponding to temperatures 10 6 K to 10 9 K, is 
considered. Exponential growth of the r-mode is found to cease at small amplitudes w 10 -4 . For 
strongly damped, low temperature models, a few modes dominate the dynamics. The behavior of 
the r-mode is complicated, but its amplitude is still no larger than about 10 -4 on average. For high 
temperature, weakly damped models the r-mode feeds energy into a sea of oscillators that achieve 
approximate equipartition. In this case the r-mode amplitude settles to a value for which the rate 
to achieve equipartition is approximately the linear instability growth rate. 

PACS numbers: 04.40.Dg, 04.30.Db, 97.10.Sj, 97.60.Jd 



In a previous paper Q] (henceforth Paper I), we pre- 
sented details on the computation of the coupling co- 
efficients and damping/driving rates for the generalized 
r-modes of a uniform density, uniformly rotating star in 
the limit of slow rotation. In this paper we explore the 
nonlinear dynamics of the oscillator network constructed 
in Paper I. We concentrate on the results that will have 
physical implications for the saturation of the r-mode in- 
stability in neutron stars. The astrophysical motivation 
for this problem was given in Paper I. 

Consider a sea of coupled oscillators obeying the am- 
plitude equations 

C A (t) - iw A C A + 1ACA = -i— ^2 K ABC°B C C , (1) 

€A BC 

where the coupling coefficients k abc the driv- 
ing/damping rates ja and the frequencies wa have been 
derived from the fully nonlinear problem by means of sec- 
ond order perturbation theory. Only the most significant 
nonlinear and driving terms have been retained. Thus, 
we assume that modes never achieve large amplitudes 
c a , and explore the dynamics in the limit of weak non- 
linearity. The validity of this assumption can be assessed 
in retrospect. 

The oscillator network for the r-mode problem has 
many distinct properties inherited from the original phys- 
ical problem: The rescaled frequencies of the oscillators 
are bounded between -1 and 1, the connectedness of the 
coupling graph is determined by the selection rules and 
the general background of the coupling coefficients pro- 
vides the landscape on which the dynamics will take 
place. 



In neutron stars energy enters the system of coupled 
oscillators at a large scale via a comparatively weakly 
connected r-mode. The energy is drained out of the sys- 
tem by viscosity at smaller scales, from modes that are 
comparatively well connected and have several couplings 
to other small scale modes. The passage of energy from 
the large scale to the damping scale is what we propose 
to study as a mechanism for saturating the r-mode insta- 
bility. 

Many examples of driven fluid systems exist in nature; 
Rayleigh-Bernard convection Q and Taylor- Couette flow 
are two particularly well-studied examples. Both dis- 
play a rich pattern-forming transition from laminar base 
flow to more turbulent regimes. The r-mode problem can 
be expected to display similar complexity. The experi- 
ments that follow are designed to indicate the types of 
behavior that may arise in different parameter regimes 
representative of those in a neutron star. 

We begin by considering the simplest coupled system 
of three modes in Section [I] Using ideas obtained from 
the three mode problem we compute possible paramet- 
ric instability thresholds that arise from the nonlinear 
coupling of the inertial modes of an incompressible 
uniformly rotating fluid. We investigate the effect of tem- 
perature on the thresholds of the inertial modes of the 
three mode system. We use the analytic results for the 
instability thresholds and oscillation periods obtained in 
Section[l]to check our integration scheme and to indicate 
which amplitude regimes are important. 

The integration of approximately 5000 complex oscil- 
lators, in which all the inertial modes with n < 30 are 
included, coupled by approximately 1.3 x 10 7 coupling 
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coefficients, for integration times of > 10 7 /2f2 requires 
equation to be implemented as efficiently as possi- 
ble. Our approach to implementing Q is detailed in 
Section HJ 

We build up our intuition about the behavior of 
our coupling network in Section IIIII by integrating the 
triplet of modes in our network that go unstable first, 
namely the triplet having the lowest parametric instabil- 
ity threshold. A second and third direct coupling to the 
r-mode are then systematically added. Any simulation 
starting from negligibly small initial conditions will dur- 
ing its initial evolution display the dynamics observed for 
these low order systems. 

In Section |V| we begin our investigation of the large 
modal system by observing how rapidly if at all equipar- 
tition of an initial large excitation in the principal r-mode 
takes place. The opposing concepts of ergodicity and 
phase coherence are introduced, and an indication of the 
nonlinear timescales is established. 

Two damping regimes exist. When dissipation is 
strong, the daughter modes that couple directly to the 
r-mode can damp out the r-mode's amplitude growth 
effectively. The daughter-daughter couplings may alter 
the dynamics but are not essential in halting the growth. 
When dissipation is weak, the modes to which the r-mode 
goes unstable do not damp its growth effectively and mul- 
tiple generations of daughter-daughter couplings are nec- 
essary to transfer the energy to the damping modes. To 
explore this case fully we need to introduce more cou- 
plings. However, it is impractical to include all the cou- 
plings since their number scales roughly like ~ n 8 . A 
criterion for identifying the important couplings is essen- 
tial to render the calculation tractable. Our method of 
dealing with the proliferation of couplings is discussed in 
Sections HU and IVTBl 

We conclude our investigation by observing the dy- 
namics of a few weakly damped systems in the high tem- 
perature regime. In this regime the transition to more 
Kolmogorov-like spectra is made and a large number of 
modes are excited. 



I. THREE MODE SYSTEMS 

The smallest non-linear system consists of three modes 
and one coupling. The undamped Hamiltonian system 
yields an analytic solution in terms of elliptic integrals 
[5j as shown in Appendix ^ The period of the am- 
plitude oscillations depends on the initial conditions of 
the system and can be accurately calculated numerically 
0. For the case of a large amplitude parent coupling 
to small amplitude daughters the period increases as the 
parent amplitude decreases. It was this dependence on 
initial amplitudes of the nonlinear oscillation frequency 
that was used to explain the timescales involved in 
the catastrophic decay of the "pumped up" state of the 
r-mode observed in hydrodynamical simulations ||. 

There exists a minimum parent amplitude below which 



no oscillations in amplitude occur. This amplitude is the 
no-damping limit of the parametric instability threshold 
discussed in the next paragraph. Above this threshold 
very long oscillation periods are realized. Accurately fol- 
lowing these oscillations, which take place over thousands 
of periods of the oscillators involved, is a challenge to our 
numerical code, and was one of the tests used to verify 
the reliability of our integration scheme. 

Study of the damped system introduces the concept of 
the parametric instability, which sets the minimum am- 
plitude the r-mode has to reach before it can excite other 
modes in the system. No energy is transfered out of the 
driven mode until the instability threshold is reached. 
The phase space of the damped three mode system was 
thoroughly explored by Dimant |9j and Wersinger et. al 
|10|. The stationary solutions and parametric instabil- 
ity threshholds are compactly re-derived in Appendix IBl 
The parametric instability threshold for an unstable par- 
ent a sharing excitation with two daughters (3 and 7 is 



7/37 7 

AW-yWpK 2 



Sw 



77+7/3 



(2) 



Here 7^, 7 7 are the damping rates and wp, iu 7 are the fre- 
quencies of the daughter modes, while Sw = w a —wp — w^ 
is the detuning. Note that the instability threshold does 
not depend on the linear growth rate of the parent and 
is subject to the correct relationship between the signs of 
the mode frequencies (some triplets never go unstable as 
shown in Appendix IBl) . 

The parametric instability thresholds for the 146,998 
direct couplings to the r-mode were calculated for a range 
of temperatures. All modes with n < 30 were included. 
The behavior of the lowest three thresholds as a function 
of temperature is shown in Figure ^ 

The graph begins with the temperature at which grav- 
itational radiation first destabilizes the r-mode of a neu- 
tron star rotating at fl = 0.37^/irGp, namely where the 
balance between the damping due to shear viscosity of 
the r-mode and the GR driving rate occurs. At this tem- 
perature all modes have strong viscous dissipation. The 
r-mode has to attain a relatively large amplitude before 
the nonlinear term dominates the linear damping of the 
daughters (equation J2J). By contrast in a system where 
the daughters are less strongly damped parametric insta- 
bility occurs more readily. As the damping is decreased 
the instability amplitudes decrease monotonically until 
the limiting value of \Sw / Ak^/w^w^\ is reached. 

The inset in the graph gives the conversion of the num- 
bers labeling a particular coupling to the eigenmodes la- 
beled according to degree, n, and order, m, of the Leg- 
endre functions that make up the eigenmodes' functional 
representation. The lowest parametric instability in the 
network we considered occurred for modes with reason- 
ably high degree, namely n = 13 and n — 14. For low 
enough temperature, the shear damping is large enough 
so that the three mode system can dissipate all the energy 
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FIG. 1: Temperature dependence of lowest 3 parametric instability thresholds of the r-mode to the inertial modes. 



with which the r-mode is driven, stopping its growth. 
However, as the damping decreases the r-mode is driven 
at a rate that cannot be halted by mere three mode cou- 
pling, and more modes become important. 



II. NUMERICAL IMPLEMENTATION OF THE 
AMPLITUDE EQUATIONS 

We implement the amplitude equations Q numerically 
by means of a Bulirsch-Stoer method with adaptive step- 
size control 0. Adaptive stepsize is a feature that is es- 
sential to any scheme adopted. The evolution begins with 
slow exponential growth described by linear perturbation 
theory, during which there is little excitation of the rest 
of the network, allowing the use of larger stepsizes, but 
proceeds to the excitation of many modes via nonlinear 
effects, which demands shorter stepsizes. A further cri- 
terion for the selection of an integrator is that it must 
minimize the number of function evaluations of the right 
hand side (RHS), since looping over a large number of 
coupling coefficients becomes very expensive. 

Another integration scheme that was explored was the 
fourth-order Runge-Kutta scheme, also with adaptive 
stepsizing 6]. However, this method displayed a larger 
systematic secular drift in the energy of a Hamiltonian 
system than the Bulirsch-Stoer method did. Symplec- 
tic integrators were also investigated, but we were not 
able to find one with adaptive stepsizing that could com- 
pete with the high order accuracy of the Bulirsch-Stoer 
method with roughly the same number of function eval- 
uations. Furthermore, as soon as dissipation is added 
the advantage of using the symplectic integrator is no 
longer clear. Semi-implicit extrapolation methods, such 
as the stifbs routine Q , allow larger stepsizes to be taken 
at the cost of the construction of the Jacobian matrix. 
Although equation Q is well suited to the construction 
of the Jacobian analytically, limitations on memory, and 



the cost of communicating large matrices in a parallel ar- 
chitecture, made this integrator viable only for a system 
with a small number of modes. Regardless of the inte- 
grator used, the majority of the CPU time will be spent 
in the evaluation of the RHS of equation Q). This is a 
major limiting factor in integrating over very long times. 

To implement the RHS of equation^ we maintain an 
ordered list of all possible nonzero couplings ^xbc ^ or 
each entry in the list we store the three modes that are 
coupled, their coupling coefficient, and the frequency de- 
tuning. The modes are identified by the one dimensional 
index j introduced in Paper I which arranges the iner- 
tial modes labeled by the three numbers (n,m,fc) con- 
secutively. A vector c contains the complex amplitudes 
indexed by the same index j. 

After initializing the derivative terms to zero (c = 0), 
we loop through all the couplings, adding the contribu- 
tion from each nonlinear three mode coupling to the RHS. 
Since the amplitude vector is indexed by the mode num- 
ber this can be done very efficiently (Note that all the 
coupling coefficients for the problem we are considering 
are real.) 

Two possible cases exist. If j3 = 7, then the rule for 
updating the amplitudes due to this coupling is 

C/3 — > Cp - i2wpKappC a C*p . (3) 

Otherwise (the general case) , when the coupling is among 
three distinct modes, the rule for updating the RHS is 

C/3 — ► Cp - i2lUpKaP 7 C a C* 

Cry !" C-y l^UJy K-QiPyC a C p . (^) 

The factor of 2 takes into account the symmetry in the 
last 2 components of the coupling. (See [llj . equation 
(4.17). Note that for our modes, the m-selection rule 
forbids the last term of equation (4.17) from occuring.) 
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Constructing the RHS in this manner lends itself to 
parallelization especially when systems with large num- 
bers of coupling coefficients are integrated. Distributing 
the coupling coefficients over the various processors, let- 
ting each processor compute the contribution due to its 
local set of coupling coefficients, and then combining the 
result by communicating only the relatively small deriva- 
tive vector results in a very efficient parallel code, with 
only a small overhead in communication costs. The com- 
putation time decreases roughly linearly with the number 
of processors added for systems with a large number of 
coupling coefficients. 

After communication between the processors has taken 
place, the linear part of the equation is added using 

+ (iw a - j a )c a . (5) 

This implementation follows each individual oscillation 
of the real and imaginary parts of a modes involved, ft 
yields very accurate results for energy conservation of the 
Hamiltonian system, 

E = ^C Q C*-4 22 Ka/3 7 K(c*C / 3C 7 )-2 ^ Ka/3/3^(c*c|) . 

(6) 

In equation © summation is implied over all nonzero 
couplings, (not mode indices as used in .1. l,j ) and 
indicates the real part of a complex number. 

For the scheme outlined above, the accurate integra- 
tion of the frequency terms iw a c a limits the timesteps 
to less than a rotation period, making it impractical 
for achieving long integration times of approximately 
I0 7 /2O. The coordinate transformation 

c a = C a e ° (7) 

removes the explicit oscillatory term. However, it does 
so at the expense of introducing a few rapidly oscillat- 
ing terms, proportional to e lSwt (see equation iJHUO m 
the nonlinear coupling terms. These terms may have fre- 
quencies up to three times those found in the original sys- 
tem, equations QJ. Thus the coordinate transformation 
itself results in no significant speed up of the integration. 
However, the couplings that play the most important role 
in the long term evolution of the problem are those near 
resonance, with detuning close to zero. Couplings with 
large detuning that oscillate very rapidly limit the step 
size of the simulation, but on average contribute very 
little to the longer term dynamics. Retaining only cou- 
plings with a detuning less than some pre-selected cutoff 
allows much larger timesteps to be achieved. For ex- 
ample retaining only couplings with Sw < 0.002 for the 
n < 30 system allows stepsizes w 400/2S7 to be taken 
initially and as more modes become active this decreases 
to about 40/2f2 . The validity of this assumption and the 
effect of choosing different cutoffs are explored in greater 
detail in Section IYl Bl 

Imposing an upper limit on 5w also decreases the 
number of couplings retained in the nonlinear interac- 
tion terms, thus speeding up the calculation. Additional 



speedup is achieved by dynamically adding couplings as 
the simulation progresses. Consider an initial condition 
in which all the modes lie below a certain noise thresh- 
old, chosen small enough that the non-linear interactions 
are negligible. Since at most a few r-modes are unsta- 
ble at first, the evolution is simply exponential growth 
of these modes. As these unstable modes grow in am- 
plitude, they drive the modes that couple to them di- 
rectly. All couplings that do not couple directly to the 
unstable modes during this early phase can be neglected. 
As an unstable mode passes its first parametric insta- 
bility threshold, the daughter modes coupling directly 
to this mode will increase in amplitude, rise above the 
noise, and start playing a significant role in the dynam- 
ics of the simulation. When the daughters rise above a 
preselected threshold their contribution to the nonlinear 
driving will become important for other modes, and all 
direct couplings to the daughters are then added to the 
existing coupling list. The step size is decreased and the 
integration continued with the updated coupling list. In 
practice, to prevent repeatedly restarting the simulation 
with new couplings, when a certain mode is tagged as be- 
coming unstable, the direct couplings to all latent modes 
with an amplitude within 50% of the amplitude of the 
unstable mode are also activated. Successively adding 
couplings in this manner allows the simulation to move 
through the initial phase of computation more rapidly. 

We used a number of checks on our integration scheme. 
Analytic results for the three mode Hamiltonian system 
are reproduced, including the sensitive dependence of the 
nonlinear period for parent amplitudes just above the 
parametric instability threshold. The parametric insta- 
bility thresholds for the damped three mode system are 
accurately observed. (These thresholds continue to be 
used as an analytic guide for larger systems although 
there is no strict mathematical proof that ensures the 
modes should become unstable at the lowest three mode 
parametric instability.) We also checked that for zero 
damping and growth rates the energy of the Hamiltonian 
system © drifts very slowly compared to the growth rate 
of the r-mode instability. 



III. INTEGRATIONS OF SMALL SYSTEMS 

As was mentioned in Section [I] and Appendix [0 the 
damped three mode system has been studied extensively 
El [HI 03 ■ However, there have been only a few attempts 
at quantifying the dynamics in the case where the steady 
state solution is unstable, for example, when the detun- 
ing is small or the daughter mode damping is slow. As a 
result it is instructive to have a closer look at the partic- 
ular three-mode coupling in our network that goes unsta- 
ble first for various damping rates, as its dynamics sets 
the initial development of any simulation. The existence 
and character of steady state solutions and the behav- 
ior of the three mode problem when the daughters can 
no longer effectively damp out the parent govern which 
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FIG. 2: Numerical solutions to the three mode problem with the lowest (7/3,77 — > 0) parametric instability threshold for a 
variety of temperatures (a) T = 5 x lO 6 ^ (Chaotic motion), (b) T = 9 x lO 6 ^ (two point periodic cycle), (c) T = 1.5 X 10 7 .K" 
(approaches one point periodic cycle) and (d) T = 2 x 10 r if (diverging solution). The actual ( 7/3, 7 7 7^ ) parametric 
instability thresholds are indicated with dashed lines. 



modes are excited subsequently. 

The evolution of the three-mode system with the low- 
est (zero damping) parametric instability threshold is 
shown in Figure [21 for various damping rates relevant to 
our problem. 

The qualitative picture in which the parent mode is 
damped by two daughter modes and reaches a steady 
state solution presented by Arras et al ^| is somewhat 
simplistic. Rather, the parametric instability threshold's 
strong dependence on frequency detuning (see equation 
implies that a coupling with small detuning will go 
unstable first. Even if this coupling is strongly damped, 
the small detuning places it in a regime (as clearly de- 
picted by Dimant in Figure 1(b) of his paper where 
no constant amplitude states exist, and all solutions vary 
with time. Furthermore, the closer to resonance the cou- 
pling is, the more complex the solutions become, dis- 
playing increasingly many periods and entering chaotic 
regimes [ljj. The first unstable coupling of the n = 3, 
m = 2 r-mode falls in this regime, so large unpredictable 
changes in its amplitude may result, perhaps causing the 
r-mode to be hard to detect conclusively in the strong 
damping regime. 

In the plots (a) to (c) of Figure [5] the daughter modes 
are sufficiently damped to saturate the growing r-mode 
via a single three-mode interaction. The maximum am- 
plitude the r-mode reaches is roughly of the same or- 
der of magnitude as the parametric instability threshold. 
As shown in Figure ^ large damping has the counter- 
intuitive effect of resulting in large amplitudes. One of 
the identifying features of these strongly damped systems 
is the saw-tooth shape of the modal evolution: exponen- 
tial growth of the parent, followed by excitation of the 
daughter modes, and rapid dissipation. This saw-tooth 
feature is also seen in evolutions with large numbers of 
modes. As the damping is decreased, the nature of the 
three mode solution changes from a system bounded in 
phase space and energy to a system that is diverging 
with ever increasing frequency oscillations, Figure [3d). 
High temperature neutron stars most likely fall into this 
regime. 

The energy for diverging three mode systems is shown 



in Figure |3Ja) for various temperatures. A period of ex- 
ponential growth increases the parent amplitude above 
the parametric instability threshold, whereupon energy 
begins to be transfered to the daughter modes. Subse- 
quently the total energy diverges at a slower rate. The 
divergence rate increases roughly linearly with decreas- 
ing daughter mode damping. High temperature simula- 
tions diverge more quickly than the more damped low 
temperature simulations. Introducing the couplings with 
the second and third lowest parametric instability thresh- 
olds decreases the growth rate of the Hamiltonian energy 
further, as is shown in Figures Bib) and[2Jc). Empiri- 
cally, the dynamics begin to differ from the unstable three 
mode evolution once the parametric instability thresh- 
olds of the additional mode couplings are exceeded. The 
instability thresholds of the couplings used to make these 
plots as well as the modes involved are shown in Figure^ 
The first and second parametric instabilities occur for 
relatively low order modes with n ~ 13, 14, 15 where the 
viscous damping cannot stop the growth entirely. The 
third instability is to modes with n ~ 23, 24, where dissi- 
pation is large enough to stop the growth, so the system 
energy levels off. 

For the most rapidly diverging, T = 3.7 x 10 7 K case 
shown in Figure EJa), the isolated three mode solution 
has wild oscillations with ever increasing frequency and 
amplitudes. Adding the second coupling Figure 0tb) 
changes the dynamics dramatically: After the second 
parametric instability is exceeded, the second coupling 
dominates, while the initially unstable three mode cou- 
pling is almost entirely suppressed. The third coupling 
Figure ^c) becomes important at an amplitude between 
the parametric instability thresholds for the second cou- 
pling and the third couplings. When the r-mode reaches 
the third parametric threshold, the dynamics change 
once again; thereafter, it appears as if the second cou- 
pling no longer suppresses the first triplet, whose daugh- 
ter modes begin to grow again. 

General three mode systems have been extensively 
study by Baesens et al [lj| and can display very com- 
plicated behaviour in various regions of parameter space. 
Adding additional oscillators further increases the phase 
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FIG. 3: Evolution of the Hamiltonian energy for temperatures in the range T = 2.1 x 10 K to T = 3.7 x 10 K in increments 
of 2 x 10 6 K. Plot (a) shows the diverging three mode system with the lowest parametric instability. Plot (b) shows the system 
used in (a) with the coupling that yields the second lowest parametric instability threshold added. Plot (c) shows the case with 
the three lowest parametric instability thresholds included. Dashed horizontal lines indicate the linear energy corresponding to 
the parametric instability thresholds. 
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FIG. 4: (Color online) Evolution of the modal amplitudes for the systems discussed in Figure[3]for the case with T = 3.7 x 10 7 K 
Dashed lines indicate the parametric instability amplitudes. The n=3, m=2 r-mode is indicated in black, other colors are chosen 
to match the color coding for the instability thresholds used in Figure 



space to be explored, and the complexity of solutions that 
can be attained. Our study of these systems is not com- 
prehensive, and only serves as one example of what can 
happen as additional degrees of freedom are added. It is, 
however, the example that is most relevant to our larger 
network. One can hope to recognize the effects of indi- 
vidual three mode couplings in the larger systems. This 
can allow one to differentiate between the effects of indi- 
vidual three mode couplings and those that result from 



the combined couplings of many modes. 

The previous examples do not take into account the 
couplings among the daughter modes which will be 
present in a larger system. As we turn to coupled sys- 
tems with numerous modes, additional factors such as 
the effective number of degrees of freedom in the phase 
space available at a certain time, and the structure of the 
coupling network, will become important. 
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IV. LARGER SYSTEMS 

During subsequent sections we explore the character- 
istics of an oscillator network that includes all of the in- 
ertial modes up to n — 30, with the gravitational driving 
fixed to the value for a star rotating at = 0.37y/irGp 
(i.e. v — 695 Hz). This corresponds to a principal 
(n = 3, m = 2 or j = 4) r-mode driving rate of 
74 = 5.9 x 10~ 7 (2O) in our dimensionless units. Including 
all couplings among these modes is impractical. There 
are a total of approximately 10 9 nonzero couplings. So we 
only account for couplings among modes with frequency 
detunings below a pre-selected maximum value. The ef- 
fects of employing different detuning criteria are explored 
in greater detail in Section IVIBI The majority of sim- 
ulations are integrated using a criterion of Sw < 0.002 
which results in the inclusion of 1306999 couplings be- 
tween the 4995 modes. We begin by integrating the 
Hamiltonian system; this checks our integration scheme, 
and also yields insights on the interaction among numer- 
ous modes. We then continue to explore our large net- 
work for various shear damping regimes in Sections I VII 



V. LARGE HAMILTONIAN SYSTEMS 

Without damping and driving, the Hamiltonian en- 
ergy of the oscillator network (equation ©) should be 
conserved for any subset of couplings. Having an in- 
tegration scheme that accurately conserves energy over 
a long period of time is particularly important because 
otherwise the slow growth and damping rates of the r- 
mode problem may be overshadowed by spurious numer- 
ical effects. Figure [S] (top) shows the fractional change 
of the Hamiltonian energy (E) for a system in which the 
r-mode is initially excited to an amplitude just above the 
second parametric instability threshold. Note that ini- 
tially E = 1.69 x 10~ 8 . The fractional energy change 
was 4 x 10~ 7 over a time period of 9 x 10 6 /2S7; in the 
same time span, exponential growth of the r-mode would 
result in AE/E ~ 10. Although the growth rates in 
a large nonlinear system may be much smaller than in 
linear theory, we feel confident that our simulations can 
accurately track AE/E > 10~ 5 on timescales ~ 10 7 /2S1 
for fully nonlinear systems including the gravitational ra- 
diation instability and viscous decay. 

A large Hamiltonian system involving three mode cou- 
plings that has been extensively studied is the Fermi- 
Pasta-Ulam problem (FPU) 0] ^(| . First performed in 
1954, this numerical study of a chain of coupled nonlinear 
oscillators yielded a number of surprises. Fermi et al. ex- 
pected to observe relaxation to an equipartitioncd state, 
when starting from an initial state in which only a large 
scale mode was excited. Instead, they observed that only 
a few modes were excited; moreover after a long enough 
time, the system appeared to return to approximately the 
initial state (recurrence). The remaining modes never 
rose out of the noise. Subsequently the FPU problem 
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FIG. 5: (Color online) Evolution of the Hamiltonian system 
including the 1306999 couplings with 8w < 0.002. An ini- 
tial excitation of 1.3 x 10 -4 is placed in the principal r-mode. 
n — m + 1 r-modes are indicated in black, all other modes are 
indicated in light grey/cyan on the lower panel. The paramet- 
ric instability thresholds are indicated with dashed lines. The 
relative energy change durning the simulation is indicated on 
top. 



has been studied extensively 01 \lM E3 • Because the 
dispersion relation of the FPU problem admits very few 
direct resonances among modes, in retrospect one might 
have expected equipartition to fail to be achieved for the 
small initial amplitudes chosen by Fermi et. al. Although 
the FPU problem does tend to equipartition for initial 
amplitudes above a certain amplitude threshold |2l|, its 
failure to do so at small amplitudes cautions us that a 
Kolmogorov spectrum |22( (of which the equipartition so- 
lution is a special case), may never arise in a system in 
which relatively little energy is fed into a mode that has 
few direct resonances to other modes. 

In the network we are considering the r-mode has very 
few near resonances. The nearest resonance, with a de- 
tuning of 1.9 x 10~ 6 , is associated with the coupling that 
has the smallest parametric instability threshold in Fig- 
ure n The small scale daughters to which the r-mode 
goes unstable first are directly coupled to many more 
modes. The dispersion relation for inertial modes, which 
implies that at a given n there are at least n different 
modal frequencies roughly uniformly distributed between 
— 2fi and 20, ensures that small scale modes may have 
many near-resonant couplings. After the initial bottle- 
neck to energy flow out of the large scale n = 3, m = 2 
r-mode, equipartition among the daughter modes is much 
more probable than in the FPU system. The lower panel 
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of Figure [5] displays the decay of an r-mode with initial 
amplitude just larger than the second parametric insta- 
bility threshold to a more "equipartitioned state" . At 
late times, the r-mode amplitude appears to approach 
the first parametric instability threshold, rather than a 
true equipartitioned state. However, we did not run our 
code long enough to verify that the r-mode amplitude 
never falls below the first instability threshold. 

A measure of the degree of equipartition of our system 
as a function of time is the spectral entropy p3j 



S(t) = -£>,-(*) In wj(t) 



(8) 



where the weights Wj are given by the fraction of linear 
energy in the jth mode: 



Wj 



(9) 



In equilibrium S(t) reaches a maximum value of S rnax — 
lii(N mo d es ) where N mo d es is the total number of modes 
coupling. In this state the system is likely to become 
ergodic and could be described by the random phase ap- 
proximation. In order to compare the evolutions of two 
different systems, the entropy can be normalized using 



s(t) 



S{t) - S„ 
S(0) - s v 



(10) 



where s « 1 when the energy is restricted to a small 
number of modes, and s = for equipartition. Figure |H| 
shows s(t) for three different initial r-mode amplitudes, 
chosen to lie between successive parametric instability 
thresholds. For larger initial energies, the decay to an 
equipartitioned state occurs rapidly. As the initial am- 
plitude is decreased, the timescale to reach equipartition 
increases. A system with driving and damping may re- 
main far from equilibrium, with dynamics dominated by 
just a few modes if the time to equipartition exceeds 
the damping and driving timescalcs. In this situation, 
the random phase approximation fails, and the evolu- 
tion of the system may be rich and varied. The pattern- 
forming transition to turbulence, which is well-studied 
in Rayleigh-Bcrnard convection, is an example of this 
regime 2] . 

For our oscillator network, ergodicity begins to fail for 
initial r-mode amplitudes between c a (0) = 0.8 x 10~ 4 
and c a (0) = 1.3 x 10~ 4 . The approach to equipartition 
for c a (0) = 0.8 x 10~ 4 is very slow. As is shown in 
Figure [7J a small subset of modes maintain apparently 
phase-coherent dynamics, reminiscent of the FPU prob- 
lem. 

An alternative mechanism for damping the r-mode 
growth is nonlinear interactions with a sea of coupling 
modes. The observed decay toward equipartition can 
be considered as an effective many- mode damping effect. 
For sufficiently large amplitudes, details of the network 
appear to be insignificant and the time to equilibration 
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FIG. 6: (Color online) Approach to equipartition of the 
Hamiltonian system. The time evolution of the spectral en- 
tropy for various initial r-mode amplitudes : The dashed line 
is c Q (0) = 5 x 10" 4 (blue), the solid line is c a (0) = 1.3 x 1(T 4 
(cyan) and the dot-dashed line is c Q (0) = 0.8 x 10 _4 (red). 



rapid. In the absence of dissipation, the energy originally 
in the r-mode will be shared rapidly among the network 
of coupled modes. In a very large system of modes, the 
"specific heat" is roughly proportional to the number of 
modes. As a result, the amount of energy needed to 
raise the effective amplitude of the coupled system is very 
large. In a system that has achieved equipartition, the 
amplitudes of all modes are (statistically) the same so the 
energy cost to raise the amplitude of any single mode - 
including the n = 3, m = 2 r-mode - by a specific amount 
is increased by the number of modes in nonlinear contact. 

As the initial r-mode amplitude and energy in the sys- 
tem are decreased, the details of the network become 
more significant, governing which modes can be excited 
and what pathways the energy transfer follows. The time 
to equipartition becomes larger as the amplitudes de- 
crease, effectively trapping the energy in a few modes. 
If the r-mode is driven by gravitational radiation start- 
ing from very small amplitudes, its energy will change 
rapidly at first compared with the timescale to equipar- 
tition. However, as its amplitude rises further, the r- 
mode will begin to share energy with other modes, so 
that its amplitude will increase much more slowly subse- 
quently. Thus, the effective rate of growth of the r-mode 
instability will vary as the mode amplitude evolves, but 
eventually settles down to a rate slower than indicated 
by linear theory as a result of nonlinear interactions. 

In Section IIIII the concept of decreasing the r-mode 
growth rate due to viscous damping of one or several 
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FIG. 7: (Color online) Modal dynamics for a Hamiltonian 
system with an initial amplitude c a (0) = 0.8 x 10~ 4 (bot- 
tom). The n — m + 1 r- modes are indicated in black. The 
distribution of the maximum modal amplitude obtained with 
respect to mode number is shown on top. 

daughter modes was introduced. In a many-mode sys- 
tem, the effect of dissipation is more complex than just a 
drain of energy. Large damping increases the parametric 
thresholds, thus limiting the pathways available for en- 
ergy flow, and effectively decreasing the "specific heat" 
or many-mode damping effect, which results in a greater 
growth rate of the r-mode. 

From an observational point of view, if an equiparti- 
tioned solution or a more general Kolmogorov solution is 
achieved, gravitational radiation would be emitted over 
a continuous spectrum from the sea of inertial modes, 
rather than a discrete spectrum at the r-mode frequency. 



energy effectively, and the r-mode can damp almost to 
the noise level before the daughter modes themselves de- 
cay. At such strong damping, only a small number of 
modes are important to the dynamics, and the behaviour 
is qualitatively very similar to the saw-tooth chaotic mo- 
tion observed in the strongly damped three mode case 
shown in Figure [21 

The integrations of the large network of modes cor- 
responding to the temperatures used to produce Fig- 
ures El^a) - |2Jd) are shown in Figure |H1 In these simu- 
lations only a small number of modes attain significant 
amplitudes, while the rest simply decay via shear vis- 
cosity. The amplitude spectrum from the results in Fig- 
ure Hfb) is compared with the equipartitioned spectrum 
obtained for zero damping with c Q (0) = 5 x 10~ 4 in Fig- 
ure |5J Both simulations start with an initial noise level 
of -\/2 x 10~ 6 (indicated by means of a dashed line in Fig- 
ure^ . The imprint of the shear damping is clearly visible 
in the damped spectrum (compare with the analytic re- 
sults in Paper I). The equipartitioned result just increases 
the energy in every mode by the same amount causing 
the spectrum to rise above the noise level, whereas there 
is distinct structure in the power spectrum for the dissi- 
pative system. These spectra indicate the two extreme 
cases that could exist in our simulations. 

Note that the the set of modes to which the r-mode in 
Figure IHl^a) goes unstable first is different from the first 
unstable triplet in the Figures Efb)- IHfd) and the triplet 
used to produce in Figure [S^a). The parametric insta- 
bility to which it does go unstable is correctly predicted 
in Figure G] The difference between the modes excited 
in Ela) as compared to those excited in Efb) and IBIc) 
indicates the influence the first unstable triplet has on 
the subsequent second generation daughter modes. Fig- 
ures|H{b) andOJc) are very similar to the pure three mode 
problems depicted in Figure □» and|3c); this is to be 
expected since the damping is strong enough to saturate 
the r-mode using the three mode coupling only. Fig- 
ure |Efd) corresponds to a diverging solution in the three 
mode problem Figure|2Id) . This case shows a larger num- 
ber of modes raising above the noise. 



VI. LARGE NON-CONSERVATIVE SYSTEMS 

A. Strong Damping Case 

The inclusion of driving and damping terms adds new 
nondynamical timescales to the problem, and changes 
the instability thresholds. An estimate of the time scale 
introduced by gravitational radiation for our model is 
l/7 4 = 1.7 x 10 6 /(2f2) . For very strong damping, the 
instability threshold of the daughter modes increases, so 
the r-mode reaches large amplitudes before the daugh- 
ters begin to grow. However, at these large amplitudes, 
the effective non-linear energy transfer rate is also rapid, 
resulting in rapid energy transfer from the daughters to 
higher order modes. These high order modes dissipate 



B. The effect of detuning 

In our simulations, we select the couplings that we ex- 
pect to play an important role in an evolution based on 
their detuning, because including all the couplings be- 
tween the modes up to n — 30 is impractical. Both the 
time to calculate all the couplings and the time to run 
the evolution would be absurdly long. 

The suitability of a maximum detuning criterion as a 
method of selecting important modes can be seen from 
equation (|B2(I . Modes with large detuning will oscillate 
very rapidly while a smaller detuning implies a slower 
variation and will on average have a greater effect on the 
modal evolution. Taken to the extreme, the selection of 
modes based on a detuning criterion results in including 
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FIG. 8: (Color online) Top panels: Time evolution of large strongly damped systems for temperatures (a) T — 5 x 10 6 K (b) 
T = 9 x W 6 K (c) T = 1.5 x 10 7 K and (d) T = 2 x 10 7 A". Bottom panels: Corresponding distribution of averaged modal 
amplitudes that rise above the initial noise level. In the top panels, the parametric instability thresholds are indicated with 
dashed lines. 
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FIG. 9: (Color online) Modal spectra for strongly damped 
system, T = 9 x 10 6 K (black) and an undamped system, 
"equipartition solution" (light/cyan). Initial noise level indi- 
cated with a dashed line. 



only modes that are in exact resonance. This assumption 
is inherent in the random phase approximation used by 
Arras et al |13j in their cascade solution. This simplifying 
assumption has been very successful in explainin g sp ectra 
observed in various systems of coupling waves [22]. The 
frequency spectrum for large scale modes is, however, 
rather sparse and all couplings to large modes have a 
finite detuning. The resonance assumption thus fails in 
the case of large scale modes as it will exclude all direct 
couplings to the driven modes. A finite detuning cutoff 
thus has to be selected. 



In systems where only a few modes are excited, such 
as the strongly damped cases, the parametric instabil- 
ity thresholds, equation l|B7|) . play an important role in 
determining the dynamics of the evolution. These thresh- 
olds are strongly dependent on the detuning, further em- 
phasizing the importance of near resonant couplings. In 
computing the parametric instability thresholds for all 
the direct couplings to the r-mode we found that the 
lowest three thresholds all had detuning < 3 x 10~ 4 . 

To explore the effect of selecting couplings by means 
of a maximum detuning we ran the coupling network 
for T — 5 x 10 7 K system with various cutoffs in 
Sw. The results are displayed in Figure 1101 for Sw < 
0.0001,0.0005,0.001,0.002. The smaller the maximum 
allowed Sw, the smaller the number of pathways available 
for the energy to leave the r-mode or subsequently excited 
daughter modes. A detuning criterion of Sw < 0.00005 
was found to be too restrictive, and the solution diverged 
after 10 6 /2f2. A detuning criterion of Sw < 0.0001 ex- 
cludes the coupling with the third lowest parametric in- 
stability to the r-mode and, although the evolution re- 
mains stable for > 10 9 /2f2, the dynamics differ from evo- 
lutions with less restrictive detuning criteria. The ab- 
sence of the strong damping due to the coupling to the 
triplet with third lowest parametric instability threshold 
allows the mode to grow to larger amplitudes before de- 
caying. Characteristics such as the average energy con- 
tained in the system and average r-mode amplitude are 
roughly the same as observed using a less restrictive de- 
tuning criterion. 

The benefit of selecting a strict detuning criterion is 
displayed in Figure ^] Not only does the number of 
couplings scale linearly with detuning, but the average 
timestep also increases as the maximum allowed Sw be- 



11 



8w = 0.0001 



8w = 0.0005 




.„-4 0.2 0.4 0.6 0.8 1 1.2 1.4 
1.5^ , , , , , , , 



namics will not be sensitive to inceases in the maximum 
Sw. Since many of the couplings are near resonance, the 
system displays a sensitive dependence on initial condi- 
tions, especially later in the evolution when higher order 
modes become excited. Small changes in the round-off 
accuracy and architecture of different computers results 
in modal evolutions that differ in detail (for example, the 
positions of sawtooth meanderings of the r-mode in Fig- 
ures El an d ED shift). Global properties such as the total 
energy of the system and time averaged modal ampli- 
tudes are, however, found to be robust. 
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FIG. 10: (Color online) Evolution of the (n = 3,m = 2) r- 
mode (bottom) and total energy (top) for a T = 5 x 10 7 K 
system with different detuning criteria. The parametric in- 
stability thresholds are indicated by horizontal dashed lines. 
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FIG. 11: (Color online) Comparison of the integration times 
achieved, average timestep sizes and number of couplings in- 
cluded for various detuning criteria. All evolutions were run 
on 16 processors of the NCSA Platinum cluster for a period 
of 24 hours. The temperature was chosen as T = 5 x 10 7 K. 



comes smaller. Note that because the runs for large de- 
tuning have not yet reached a regime in which many 
modes become excited, the average stepsizes shown in 
the Figure ITT1 for Sw — 0.001, 0.002 are overestimates, 
and can decrease by up to a factor of 100 as more modes 
start interacting. All subsequent runs reported here are 
done for a maximum detuning Sw — 0.002. Although 
the detailed modal evolutions may differ from that of 
the complete system, we believe that the averaged dy- 



C. Weakly damped evolutions 

As the temperature increases, the shear damping rates 
decrease, and the amplitude decay from the initial noise 
spectrum slows. The energy from the unstable r-mode 
has more time to drain into the other modes before they 
are damped, and therefore numerous modes become im- 
portant dynamically. Figure IT2"1 displays an intermediate 
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FIG. 12: Evolution of the T = 5 x 10 7 K system with Sw < 
0.002. The parametric instability thresholds are indicated by 
horizontal dashed lines 

case, T = 5 x 10 7 K, for which the ratios of the decay rates 
of the daughter modes (414 and 538) involved in the first 
parametric instability to the growth rate of the r-mode 
are 0.14 and 0.18; for the second parametric instability 
(daughter modes 494 and 592) the ratios are 0.16 and 
0.18. Crudely, we would then expect to need the r-mode 
to couple to about 10 modes with similar damping rates 
to curb its growth directly. Instead, energy is rapidly dis- 
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tributed among the daughter modes, which excite many 
more modes in turn. If one makes a plot like Figure 
for the simulation shown in Figure El then it looks very 
similar to the strongly damped case except that all the 
amplitudes are displaced upward, so that on average they 
lie at the initial noise level. Approximately half of the 
modes in the system rise above the noise, while only the 
most strongly damped modes decay below the noise level. 
This causes the energy to be shared among the accessible 
modes. In comparison with a more weakly damped case, 
T = 5 x 10 8 K displayed in Figure El larger mode am- 
plitudes are attained in certain modes in Figure IT21 even 
though the energies of the two systems are comparable. 
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FIG. 13: Evolution of the T = 5 x 10 8 K system with 5w < 
0.002. The parametric instability thresholds are indicated 
by horizontal dashed lines. Modal amplitudes are indicated 
below and the total Hamiltonian energy top. 

Figure^Jrepresents the most weakly damped case con- 
sidered during this investigation. In a realistic neutron 
star, it lies in the temperature range where the effect 
of shear viscosity becomes small, and other damping ef- 
fects such as bulk viscosity begin to dominate. In this 
slowly damped case, the parametric instability thresh- 
olds are very close to their limiting zero damping values. 
For the first parametric instability the damping rates of 
the daughters (414 and 538) are 0.0014 and 0.0018 times 
the r-mode growth rate. It would take approximately 
1000 couplings to similarly damped modes of compara- 
ble amplitude to curb the r-mode growth. Our network 
of 4995 oscillators includes enough modes with sufficient 
damping to satisfy the criteria for a stable Kolmogorov 



spectrum |22j, namely that the net damping of all the 
modes is positive. In the weakly damped regime all the 
modes rise above the noise. Energy sharing among the 
daughters is rapid, with few individual daughter modes 
rising substantially above the average amplitude. The 
r-mode, though, maintains an amplitude well above the 
ever-increasing amplitude level of the daughter modes. 
Energy tends to be "trapped" in the r-mode because it 
has few direct nearly resonant couplings. With the excep- 
tion of the unstable r-mode, the modal spectrum closely 
resembles the one observed for the Hamiltonian case, Fig- 
ure [5J with a slightly larger variance in the modal ampli- 
tudes about their mean. By the end of the simulation no 
trace of the characteristic shear damping structure could 
be seen. 

In contrast to the more strongly damped case displayed 
in Figure 1121 where the total energy peaks and then 
drops, the total energy in this system has not yet reached 
a maximum turning point. The modal amplitudes of Fig- 
ure El are expected to continue to bounce around in the 
same amplitude range as seen in the simulation so far and 
the energy is not expected to resume monotonic growth. 
This expectation is supported by simulations with a much 
more restrictive detuning criterion, Sw < 0.0001 evolved 
for 2 x 10 8 /2£1. The outcome of the evolution in Fig- 
ure !13l is less certain. The r-mode might remain at nearly 
the same amplitude ~ 10 -4 , since any increase in ampli- 
tude is severely penalized by an increase in the rate to 
equipartition, resulting in energy flow out of the r-modc. 
As can be read off Figure the time to achieve equipar- 
tition for an r-mode amplitude of 1.3 x 10~ 4 is 10 6 /2S1 
which nicely balances the characteristic driving time due 
to GR of I/74 ~ 1.9 x 10 6 /2Q . If it were to remain at 
at this fixed amplitude, the r-mode would pump energy 
into the sea of inertial modes until they reach sufficiently 
high amplitudes that the cumulative damping of the sea 
of modes halts the growth in Hamiltonian energy. At 
this stage, the modal amplitude distribution may change 
slightly to reflect the stronger damping at the higher or- 
der modes. One of the problems in obtaining this state 
numerically is that as the energy increases, and the noise 
level is raised, the nonlinear timescales of the modal os- 
cillations decrease, which cuts the timestep and there- 
fore increases integration times. The current simulation 
shown in Figure H31 took 20 days on 32 processors of the 
NCSA Platinum cluster. 



D. Rapid Rotation 

Up until this point the dynamics of a star rotating 
at £1 = 0.37^/irGp, e = 0.5 have been considered. The 
maximum rotation rate that can be realistically han- 
dled by our perturbative model of Maclaurin spheroids is 
fl = O.Qly/irGp or e = 0.81267 where the system displays 
a bifurcation to the Jacobi sequence [24| . To demonstate 
the most extreme case of rapid rotation we consider the 
damping for a star e = 0.82; this will serve as an upper 
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bound for the maximum possible driving rate that can 
be represented in our model. The temperature we use to 
determine the shear damping is T = 5 x 10 7 K (i.e. the 
intermediate damping case of Figure lT2T) . Since the time 
variable is 2£lt, damping rates become "less effective" for 
a more rapidly rotating star. The numerical values for 
the shear damping are effectively decreased by a factor 
of 0.6; the nonlinear couplings are the same for either 
in our units. In these units the driving rate of the most 
unstable r-mode is 74 = 7.49 x 10~ 6 (2f2), as opposed to 
74 = 5.92 x 10 _7 (2f2) considered previously. This implies 
that the rate to achieve equipartition will balance the in- 
stability driving rate at an r-mode amplitude ~ 5x 1CP 4 
(see Figure EJ. For such a rapidly rotating star, gravita- 
tional radiation strongly damps certain low order modes 
and dominates viscosity in the low order region. How- 
ever, the number of modes, affected is so small that it 
does not influence the dynamics at all. The evolution of 
the modal amplitudes is shown in Figure 1141 An inter- 
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FIG. 14: Evolution of the T — 5 x 10 7 K, e = 0.82 system 
with Sw < 0.002. The parametric instability thresholds are 
indicated by horizontal dashed lines. Modal amplitudes are 
indicated below and the total Hamiltonian energy in the top 
panel. 

esting feature is that the r-mode grows so rapidly via its 
linear instabilty that it passes its first parametric insta- 
bility threshhold without any permanent slowdown. This 
is because the instability grows faster than the nonlinear 
period associated with the first unstable coupling. (A 
similar dynamical effect is expected with an artificially 
ramped up radiation reaction force. If the mode grows 
too rapidly it bypasses possible instability thresholds and 
does not excite the daughter modes.) 



VII. DISCUSSION 

We have studied the effects of second order nonlinear 
couplings among inertial modes of an incompressible star 
in order to understand how they may truncate the CFS 
instability of the n = 3, m = 2 r-mode. Although the 
system we study is not a realistic version of a neutron 
star, it illustrates certain general principles that apply to 
the nonlinear evolution of the instability. 

As in our earlier studies ^30I0j we have adopted 
the philosophy that nonlinear effects can slow or stop 
the growth of the amplitude of the inertial modes of a 
neutron star so that large amplitudes are never attained. 
Thus, we can employ an expansion of the fully nonlin- 
ear hydro dynamical equations of motion, keeping only 
the lowest order nonlinear terms . This philosophy can 
be justified in retrospect, since as shown in this paper 
nonlinear effects do indeed limit the dynamics to small 
amplitudes. 

We may idealize further by restricting attention only 
to the inertial mode sector of perturbations of a (slowly) 
rotating neutron star. One reason we can do this is 
that these modes have frequencies that lie between ±217, 
which are generally far below the frequencies of the p 
and / modes of the star, which have frequencies > \JirGp 
and, for small mode wavelengths, ^> y/irGp. Thus, we do 
not expect any instability involving the r-mode to excite 
the higher frequency modes efficiently. Moreover, non- 
linear effects are magnified by near resonances even at 
small modal amplitudes. Such resonances are likeliest in 
the rather dense inertial mode sector. Indeed, one of the 
general principles governing the importance of nonlincar- 
ity in weakly nonlinear perturbation theory of the kind 
we present is that modes with small relative detunings 
dominate the dynamics. 

In this paper, we have presented a sequence of cal- 
culations of increasing complexity in order to illustrate 
the kinds of nonlinear effects that might arise. We be- 
gan by treating the simplest system: three nonlinear ly 
interacting modes. Here, a key concept is parametric in- 
stability: when the amplitude of the unstable "parent" 
mode exceeds a well-defined threshold (Equation @), it 
will drain energy into the pair of "daughters" to which it 
couples. If the daughters dissipate energy very strongly, 
then the amplitude of the unstable parent generally can- 
not grow more than a factor of a few beyond its paramet- 
ric instability threshold. However, even in such strongly 
damped cases, the dynamics can be quite complex, as is 
shown in Figurcs[2Ia)-[3Ic). Characteristically, the ampli- 
tude of the parent exhibits a sawtooth evolution and the 
system as a whole shows multiply periodic temporal be- 
havior. As the dissipation drops, though, the three mode 
system evolves toward a single period, divergent solution, 
in which the amplitude of the unstable mode increases 
secularly. In a realistic system, additional modes would 
be excited at such low dissipation rates. Indeed, adding 
the triplet interactions with the second and third lowest 
parametric instability thresholds slowed and even halted 
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the rate of growth of the instability (Figure |3J) . 

In order to understand what happens in a system in- 
volving a multitude of modes, we began by studying a 
Hamiltonian system (i.e. no unstable or damped normal 
modes). For this system, we could track the evolution 
starting from given initial r-mode amplitudes and noise 
level. We found that in general the sea of inertial modes 
tends toward equipartition, but the behavior of the r- 
mode and even a few other lower order modes is richer. 
For example, starting with an r-mode amplitude above its 
second parametric instability threshold leads to equipar- 
tition in the sea of modes, but the r-mode itself appears 
to asymptote toward its first parametric instability am- 
plitude. 

We also studied the approach to equipartition as a 
function of the initial r-mode amplitude at fixed initial 
noise level. To measure the rate of evolution toward 
equipartition we introduced the time-dependent entropy 
(equation (JSJ) ) . In general, the approach to equipartition 
is faster for larger initial r-mode amplitudes. At smaller 
amplitudes, the approach to equipartition can become 
very slow (see Figure [S] ) and may not occur at all. At 
small amplitudes the r-mode, as well as a few other low 
order modes, continue to exhibit phase coherent oscil- 
lations up to the end of our simulation at t ~ 10 7 /251. 
This maintenance of phase coherence signals a failure of 
the system to become ergodic, and also warns against as- 
suming the random phase approximation too casually in 
trying to understand nonlinear truncation of the insta- 
bility. 

Finally, we turned to multimode non-Hamiltonian sys- 
tems with growth due to the CFS instability, and damp- 
ing due to (shear) viscous dissipation. There are a num- 
ber of competing elements that determine how such a 
complicated system evolves. The two simplest are the 
growth or decay rates of the modes dictated by linear 
perturbation theory, but nonlinear interactions introduce 
several other important timescales. One was already ev- 
ident in our study of the three mode system: in the 
absence of dissipation, the modal amplitudes evolve pe- 
riodically as long as the total energy is not too large. 
A second was indicated by our study of the multimode 
Hamiltonian system: the rate of approach to equiparti- 
tion represents a kind of multi-mode dissipation rate for 
the phase coherent low lying modes. Here we interpret 
dissipation broadly as the loss of energy from large scale 
coherent motion to small scale incoherent excitations. 

For the full multimode system, we found that at strong 
damping only a few modes matter to the dynamics. The 
saw-tooth, chaotic behavior seen in the three mode sys- 
tem is also evident in strongly damped multimode sys- 
tems. Many of the modes in the sea simply decay via 
viscous damping, never having achieved substantial exci- 
tation. 

At very weak damping, though, numerous modes be- 
come involved in the evolution, which comes to resem- 
ble more closely what we observed in the approach to 
equipartition for a Hamiltonian system. Modal spectra 



resemble that of the equipartition, Hamiltonian system, 
except with larger variances of some amplitudes about 
the mean. In the weak damping regime, the total energy 
of the system continued to rise throughout our integra- 
tion time, but the r-mode amplitude appeared to level off. 
At its apparent asymptotic value, the rate of approach 
to equipartition - as indicated by Fig. 6 for our Hamil- 
tonian system - would be just comparable to the linear 
growth rate of the CFS instability. Thus, in this regime, 
the truncation of the r-mode amplitude itself is due to the 
"collective dissipation" in which energy disappears into 
a sea of coupling modes, a completely nonlinear effect. 
We expect that if we were to continue the calculations 
further in time, the energy growth would continue un- 
til the amplitudes of small scale, more highly damped 
modes became large enough for them to dissipate energy 
as quickly as the instability causes it to grow. 

The behavior outlined above for the weakly damped 
system also depends on the relative rates of the CFS 
instability, approach to equipartion, and viscous dissipa- 
tion. Wc examined this by increasing the stellar rotation 
rate, thus increasing the gravitational radiation rate. For 
this case, we did not observe any leveling of the r-mode 
amplitude, although we also did not run the simulation 
long enough (for practical reasons) for the amplitude to 
reach the value at which the growth rate of the CFS 
instability would be balanced by the rate of approach 
to equipartition. Although this does not indicate any 
inconsistency in the picture we have presented for how 
the instability is limited by nonlinear effects, it also cau- 
tions us that any calculation that enhances the instability 
growth rate could overshoot the amplitude of the r-mode 
at which it would level off in reality. 

Finally, we note that the saturation amplitudes found 
here are substantially lower than those estimated in the 
cascade picture in Arras et al. but are larger than 

but perhaps comparable to those relevant to the spin 
evolution of a neutron star estimated by Wagoner |25| . 
Our study of nonlinear effects on the development of the 
r-mode instability also suggests that the temporal evolu- 
tion of the modal amplitude is not smooth on timescales 
~ millions of stellar rotation periods. The erratic behav- 
ior seen in our simulations could, if realized in nature, 
complicate searches for gravitational radiation emission 
that depend on simple, phase coherent oscillations of the 
star. 

Because of practical limitations our modal expansion is 
truncated at n = 30. Within this truncation the r-mode 
is found to go unstable first to modes with 13 < n < 15. 
It is possible that some of the dynamics observed in this 
paper are dependent on our modal truncation. The most 
dramatic effect would result if the minimum parametric 
instablity thrcshhold were lowered. In this case the r- 
mode may go unstable to a very high order mode and 
dynamics such as equipartition rates may change con- 
siderably. It should be noted that this would imply that 
the maximum amplitude reached is smaller than or equal 
to those observed in this paper. Should the mode ever 
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grow to amplitudes comparable to those observed here, 
all mechanisms explored in the current truncation will 
still limit its growth. Our truncated expansion thus sets 
an upper bound for outcomes with larger n. The fact 
that no lower parametric threshholds were found in the 
region with 30 < n < 15 could, however, be an indication 
that we have gone far enough. Increasing the number of 
modes will have a secondary effect of increasing the "spe- 
cific heat" or the system since more degrees of freedom 
have been introduced. This may be particularly rele- 
vant to the evolution in Figure^] causing the noise level 
to rise more slowly. In a large enough system, the flat 
spectrum observed in Figure ^] should have a tail that 
is grounded below the noise level, indicating a damping 
dominated region. 

For more strongly damped regimes the modal depen- 
dence of the damping will play a significant role in deter- 
mining the accessible phase space. Sources of damping 
that are on average weak, but have a few very strongly 
damped modes will have very little effect on r-mode dy- 
namics; they will only remove a few degrees of freedom 
from a very large phase space. On the other hand damp- 
ing sources that have a fixed scaling, with large regions of 
modes damped roughly equally, will have a larger effect 
in restricting the phase space. Thus it is not the max- 
imum damping rate that is important to the dynamics, 
but rather the distribution of the most weakly damped 
modes. 



For realistic neutron stars the eigenfunction structure 
will change somewhat and this will affect the couplings to 
some extent. Changes in the eigenfrequencies could have 
a more important effect, since subtle changes could shift 
various daughter modes into or out of close resonance. 
Such changes could strongly influence the parametric in- 
stability threshholds. 

The modal frequency structure, namely that for each 
n and m = • • • n — 1 there exist n — m frequencies dis- 
tributed between — 2f2 and 2f2, should persist for more 
realistic models. It is this structure that ensures the near 
resonances among the daughter modes which in turn re- 
sults in an equipartitioned solution. The rapid increase 
in equipartition rate with increasing amplitude is also 
likely to be the result of the topology of the frequency 
spectrum and should persist in other models. 

Thanks to Larry Kidder for helping with paralleliza- 
tion of the code and numerous other technical details. 
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APPENDIX A: THREE MODE HAMILTONIAN 
SYSTEM 

Consider the amplitude equations without driving or 
damping: 

c a = iw a c a — 2iw a ncpc 1 
dp = iwpep — 2iwfjKc a c* 
dy = iWyC-y — 2iw 1 KC a dp (AI) 

Using amplitude and phase variables Cj = \/wiA i e l ^' i , 
4> = 4>a ~ 4>p — <?V 5lo — w a — wp — the system can 
be rewritten as: 

Ai = -Anti y/ WaWpw^AaApA^ sin 4> 

■ . 2nw a WRW^ cos (/>,.. . . . . x 

4> = Scu P ] == (—A a A^ - A Q Af3 + A 1 Ap) 

^/w a wpw 1 A a ApA 1 

(A2) 

where ep — e 7 = — 1 and e a — 1. The Hamiltonian for 
the three mode system with one coupling is 

H = L — 4k cos c a cpc 7 (A3) 

where L — \c a \ 2 + \cp\ 2 + |c 7 | 2 denotes the linear energy. 

Define w 3 — w a wpw 7 and Sw = Sw/w. Note that di- 
viding the amplitude equations (|A2|) by each other yields 
A a — A® a = ei(Ai — A®) , where the superscript indi- 
cates initial value at time t = 0. The difference / = L — L° 
between the initial linear energy, L°, and the linear en- 
ergy at time t can be used to rewrite the three variables 
Ai in terms of the single variable I and the three constants 
dependent on the initial amplitudes A®: 

Ai = ^l + A» (A4) 
ow 

The nonlinear term in the Hamiltonian can now be 
rewritten as follows 

where 

F(l) ={l + SwA° a )(l - SwA°p)(l - 5wA°) (A6) 

The evolution equation for I is 

I = — AnSwcaCpc^ sin0 (A7) 

or, using the Hamiltonian to rewrite sin <j> in terms of Cj 
and thus /, 

/ = T5w\J (AnCaCpCf) 2 - (L — H) 2 

= TSwsj (4k)2 (£) 3 i f, (0 - (I - H nl )* (A8) 

where H n i = H — L°. 
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Equation (|A8|) yields an elliptic integral for the time 
dependence of I that can be conveniently expressed in 
terms of Carlson's Rf function The Rf function is 
defined by 



(A9) 



where G is a cubic polynomial: 



G(t) ={t + a x ){t + a 2 ){t + a 3 ) 

=t 3 + at 2 + bt + c (A10) 

The quantities Ui(x) can be written in terms of the roots 
or coefficients of the polynomial: 

U? =U 2 + oj 



2y/cG(x) + 2c + bx 



Specializing to our problem, let 

G(l) = F(l)-Y(l-H nl f 

where 



Y = J_f6w\ 



(4k) 2 V w J 



(All) 



(A12) 



(A13) 



Thus 



G = I 3 + l 2 [Sw(A a a -A%- A%) - Y] 

+ I [2H nl Y + 5w 2 (A°A° 3 - A° a A° g - A° a A°)] 

+ 5w 3 A° Q A°A°-YH 2 l (A14) 

The time dependence can be found by integrating the 
equation l|A8(l from t = to t to yield 



1 



2VY 



-dt 



R F (UlUlUt) 



(A15) 



where Uf — Uq + U. Equation l|A15|l gives t a function 
of I, and can be inverted in terms of the Jacobian elliptic 
functions. 

The phase dependence of <fi can be found using 



or 



cos 6 = —VY 



G(0 
F(l) 

{H nl - Q 

VW) 



(A16) 



(A17) 



1. Period of Oscillation 

Theoretically the period of oscillation can be obtained 
directly from equation (1 A 15(1 . However numerical calcu- 
lation of the solution turns out to be tricky for a large 
range of energies and initial values. Scalings of the vari- 
ables and parameters turns out to be important to ensure 
accurate results in some cases. 

Introduce the parameter 



X 



~K W a Wf3W 7 

5w 2 



(A18) 



The quantity l/X is the limit of the parametric insta- 
bility threshold in the absence of damping. Then equa- 
tkmHASH can be rewritten 



l H = ±VSwXH 
Sw c 7 (0) 2 



-l H + 



Sw c a (0) 2 



H 



H 



Sw c^(0) s 



H 



W/3 



Ih + 



H 



Sw 

JTx 



Ih + 



nl 

IT 



1/2 



(A19) 



where Ih = (L — L)/H. 

Note that the Hamiltonian rescales the detuning, so 
the smaller the relevant energies the more important the 
detuning becomes. Expanding the polynomial in I in the 
square brackets in equation 1A19I) gives 



Sw 



1 H 



H 



- I'k + In— ( D a - D g - D y - — 



X 



+ ( D a Dg + D a D 1 - D g D 7 - 2 ' ' 



XSw J 



Sw 3 

"h 3 



SvH 2 i 
H 3 X 



(A20) 



where Di — Ci(0) 2 /m; 

The roots of a cubic polynomial of the form x 3 ■ 



bx 



can be computed using the intermediaries 



36 



and R = 



2a 3 - 9ab + 27c 



"9 54 

If R 2 < Q 3 there are three real roots, which is the 
only case that will be considered here. Setting 8 = 
arccos(i? / \/Q 3 ~) gives the roots as 



-2vQcos( 



9 + 2nk s 



-, k = -1,0,1 (A21) 



If the roots of the polynomial in Ih are ordered as l\ < 
I2 < h , and assuming that the / = initial conditions lie 
between roots h and I2, the period of the oscillation P is 



P 



1 



-dl 



H 



H 



-.K{k) 



y/\SwXH\ Vh^h 
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where k 2 = (I2 — h)/(h — h) and K(k) is the complete 
elliptic integral of the 1st kind. 

Note that for initial conditions with a large parent and 
two small daughter amplitudes, k gets very close to one 
(since £3 — * I2) and the integral diverges logarithmically. 



APPENDIX B: DAMPED THREE MODE 
SYSTEM 

This appendix summarizes the main properties of the 
damped three-mode system relevant to our problem. The 
various ways of rewriting the equations provide different 
possible approaches to the numerical integration of the 
oscillator net. The scheme presented by equation llB2|) 
has turned out to be particularly fruitful, especially if 
only modes of small detuning are kept. Consider the 
amplitude equations. 

Ca = iw a C a + iGRCa ~ 2iw Q KC /3 C 7 

c'p = iwpcp - jfjcp - 2iwpKc a c* 

dy = iw 1 c~ 1 — 7 7 c 7 — 2iw 7 Kc a Cp (Bl) 

Explicitly following the rapid oscilations of the individual 
oscillators can be avoided by changing variables to Cj — 
Cje~ 1Wjt , in which case the equations become 



Cp = 

C-y 



-iSwt 



- l0 Cf3 - 2i Wl3 KC a Cle lSwt 



-7-yC-y — 2iw 1 KC a C%e 



:Swt 



Using amplitude and phase variables Cj 

<t> = 4>a — <fii3 — (f>-y, Suj — w a — w/3 — w-y gives 

\c a \ =lGR\c a \ - 2w a K\cj 3 \\c 1 \ sin(0) 
M = -7/3M + 2w/3K\c a I |c 7 1 sin(^) 



<k e 



(B2) 

i<f>i 



-7 7 |c 7 | + 2w 7 k|c q ||c / 3| sin(</>) 



4> = Suj — 2k cos(0) 



-OJ0- 



I Ca 1 1 C*v 



^ [ColM 

|c/3 1 7 |C 7 | 

from which the stationary solution 

1/3 +77 _ 1GR 



t\\Cf*\ \ 



C a \ J 



(B3) 



tan© =- 



Sir 

lGRl/3 
4w a Wf3K 2 

7 7 7/3 



tan- 



tan 



(B4) 



can be obtained. 



1. Parametric instability 

If we treat the parent or driven mode as a known func- 
tion of time, then in the limit of slow driving the parent 
can be considered constant |l2j. Let d = Sie lSwt l 2 so 
that equations l|B2(l become 



C a — 
Sp = 

s* = 



constant 

- 7/3)8/3 ~ 2iwpnC a S* 
iSuj 

2iw 1 KC*S f3 + (— - 7 7 )5* 



(B5) 



Let S = [Sp S*}. Then equations iB5l> 
AS, where A is a matrix with trace Ta = —(7/3 
and determinant 



becomes S = 

7y) 



d-A = 7/37 7 + (Sw/2) 2 - Awpw 1 K 2 \C a \ 2 + (7 7 



jp)iSw/2 
(B6) 



The eigenvalues of A are A = \ (Ta ± \/T\ — 4d2) and 
the solutions are stable iff 3?(A) < where 3? indicates 
the real part of a complex number. The solution becomes 
unstable when Sft(A) = or T% = [U(y/TZ - M A )} 2 if 
Ta < 0. If Ta > the solution is always unstable. 

If a complex number u + iv = [x + iy) 2 , then u = 
x 2 — y 2 , v = 2xy or 4x — v 2 = 4x 2 u. In our problem 
u = T\ — 47^7 7 — Sw 2 + l6wpw.yK 2 \C a \ 2 , u = — 2(7 7 — 
~/p)5w and a; = at the instability point. As a result 
4Tl(-4 1 p ll -Sw 2 + 16w f3 w jK 2 \C a \ 2 ) = -4( 77 -7^) 2 ^ 2 
and the threshhold value at which the daughters or 
damped modes become unstable, known as the paramet- 
ric threshold is 



\c a \ 2 = 



77 + 7/3 



(B7) 



The signs of the frequencies are important. If in the 
above example w-yivp < 0, no parametric instability takes 
place and the zero solution remains stable for all ampli- 
tudes C' a . Suppose now that c@ is driven and the other 
two modes are damped for the same coupling. The in- 
stability criterion becomes 



M 5 



la 7 7 



Sw 



ll+la 



(B8) 



and occurs only if w 7 w a < 
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